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Abstract 


The conservative form of the unsteady Navier- 
Stokes equations In terms of vorticity and stream 
function in generalized curvilinear coordinates are 
used to analyze the flow structure of steady 
separation and unsteady flow with massive 
separation. The numerical method solves the 
discretized equations using an ADI-BGE method. The 
method _ls applied to a symmetric 12 percent thick 
Joukowski airfoil. A conformal clustered grid Is 
generated; several 1-D stretching transformations 
are used to obtain a grid that attempts to resolve 
many of the multiple scales of the unsteady flow 
with massive separation, while maintaining the 
transformation metrics to be smooth and continuous 
in the entire flow field. Detailed numerical 
results are obtained for three flow configurations 
(1) Re = 1000, a = 5°, (il) Re - 1000, a - 15°, 
(ill) Re - 10,000, a - 5°. Ho artificial 
dissipation was added; however, lack of a fine grid 
in the normal direction ha3 presently led to 
results which are considered qualitative, 
especially for case (Hi). 


1 . Introduction 

The flow over streamlined lifting airfoils has 
been a subject of considerable interest to fluid 
dynamicists and, to date, significant progress has 
been made towards the design of airfoils, wings, 
etc., by drawing together resources from 
experimental, numerical, analytical and empirical 
studies. The detailed flow structure of airfoils 
and wings near maximum lift In low-to-high Reynolds 
number (Re) flows 3till remains unresolved. The 
increasing interest in these flows stems from the 
desire for better control in civilian aircraft, to 
high maneuvering capability in high-performance 
military aircraft. The improved performance can be 
realized from the potential of Increasing maximum 
lift and simultaneously reducing drag under this 
condition. For some combination of flow 
parameters, the flow field around an airfoil 
experiences significant separation which 
deteriorates it3 performance and lead3 to the stall 
phenomenon. The nature of the stall may be 
characterized by the various phenomena such as 
separation, unsteadiness, transition and 
turbulence. The present study is directed towards 


accurately simulating thi3 flow field and providing 
further insight for this class of flows. Other 
important fluid dynamics applications involving 
unsteady flows Include blade rows in 
turbomachinery, marine propellers, helicopter rotor 
blades, and bluff bodies such as buildings, towers, 
underwater cables, etc., in cross flows. For this 
class of bluff-body flows, understanding the vortex 
shedding characteristics is very significant. The 
simulation technique presented here can also 
provide guidelines for analyzing some of these flow 
fields. 

For two-dimensional. Incompressible, unsteady, 
viscous, separated flow3 at moderately high 
Reynolds number, two viable approaches are 
available. (i) The first method Involves an 
inviscid-viscous strong- Interact ion analysis based 
on local lzed-f low regions, whereas (ii) the second 
method consists of using, in the entire region of 
Interest, a single set of equations which have the 
necessary mutual dependence between the invi3cld 
and viscous flow built into them. For massive 
separated viscous unsteady flows of interest in 
this study, the second approach is preferred. In 
this latter approach, the complete unsteady Navier- 
Stokes equations are preferred over other reduced 
forms of Navier-Stokes equations discussed by 
K. Ghia, Osswald and U. Ghia [1]. 

The prediction of two-dimensional flow past 
various airfoils has been comprehensively reviewed 
by Cebeci, Stewartson and Whitelaw [2]. Noteworthy 
papers on both the computational approaches 
mentioned earlier, as well as experimental and 
analytical methods, have been extensively reviewed 
by these authors and a3 such, no attempt is made 
here to review the past work. Instead, only some 
relevant studies using the Navier-Stokes equations 
approach and a couple of very recent experimental 
and analytical studies are cited here. Mehta and 
Lavan [3] studied the incompressible flow past a 9 
percent thick Joukowski airfoil at Re - 1000, and 
angle of incidence, a - 15°, using the Navier- 
Stokes equations, and provided accurate results for 
the stall characteristics. Care was exercised in 
determining the far-field boundary condition which 
was placed at a finite distance from the airfoil. 
Their numerical simulation employed an 0-type grid; 
the resolution of an 0-grid generally degrades in 
the wake region. The dynamic stall of an 
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oscillating airfoil was also studied by Mehta [4] 
using the Navier-Stokes equations for the N AC A 0012 
airfoil with Re up to 10,000 and the reduced 
frequency of oscillation k » 0.5 based on the half- 
chord. The formation of the leading-edge vortex 
was carefully simulated using a fourth-order 
accurate method, again, with the 0-grld. Wu and 
his coworkers [5] have contributed significantly in 
developing integro-diff erential and Integral 
methods for the solution of the Navier-Stokes 
equations for incompressible and compressible, as 
well as laminar and turbulent, flows past Joukowski 
airfoils at angle of attack. In some of these 
studies, flow past a 9 percent thick symmetric 
Joukowski airfoil at Re - 1000 and a - 15° was 
computed. In order to improve the efficiency of 
their solution, Wu et al. [6] suggested a zonal 
procedure with the unique feature of using the 
complete unsteady Navier-Stokes equations only in . 
the massive separated detached zone. For the same 
Joukowski airfoil configuration, they showed 
improved computational efficiency. Hodge et 
al. [7] solved the unsteady Navier-Stokes equations 
using O-grlds and predicted flows around NACA 

airfoils at stall Re of 4x10** and 2«10^ and 
compared their results qualitatively with 
experimental data. 

Mueller [8] experimentally studied the 
Lissaman 7769 and Miley M06~13“128 airfoils at 
Re < 300,000 which is lower than their design value 
of Re - 600,000. A hysteresis loop was observed 
which depended on the relative locations of laminar 
separation and transition. This data may become 
valuable for validating numerical analysis for 
these airfoils. Smith [9] has reviewed the 
theoretical aspects of steady and unsteady laminar 
separation. The Importance of a rational 
theoretical analysis in building a basis for the 
interacting-analysls approach and providing 
guidelines in the resolution of multiple scales by 
properly developing the grid for the Navler-Stoke3 
approach was clearly emphasized. The two- 
dimensional unsteady separation was linked closely 
with instabilities in the boundary layer and the 
separating shear layer. It is the belief of the 
present authors that coherent experimental, 
theoretical and numerical studies, with proper 
interaction amongst them, are needed to gain a more 
complete understanding of unsteady separated flows. 

The primary objective of the present work i3 
to study the two-dimensional unsteady separated 
flow past an airfoil at moderate Re near the stall 
conditions and accurately determine the flow 
characteristics in the massive separated region. 

To achieve this goal, the analysis of K. Ghia, 
Osswald and U. Ghia [1] developed earlier for 
Internal flow was used as the starting basis for 
analyzing the flow past symmetric Joukowski 
airfoils at high angle of attack. 


2. Governing Equations 

In terms of the vorticity vector u and the 

velocity vector V, the unsteady, incompressible 
Navier-Stokes equations consist of a temporally 
parabolic vorticity transport equation 

-§£ + (V-7)m = (u*V)V - — (V x 7 x m) (1) 


together with the kinematic definition for 
vorticity 

Z » 7 x V ; (2) 

Here, Re is a Reynolds number defined as 
Re - U^/v - U^C/v 

where, for the airfoil problem, the reference speed 
U R is U a , the reference length t, R is the chord c 

and v is the constant kinematic viscosity of the 
Incompressible fluid. Equations (1) and (2) have 
been nondimens ionalized using c as the 
characteristic unit of length, U m as the 

characteristic unit of 3peed, c/U^ as the 
characteristic unit of time, and U^/c as the 
characteristic unit of vorticity. 

For 2-D or axlsymmetrlc flow, the local 

velocity vector V can be related to a 3tream 
function ij; as 

v . 7* x ;3 o) 

where e^ is the fundamental contravar iant base 
vector. For consistency with Eqs. (1) and (2), Eq. 
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(3) is nondlmensionallzed using U^c as the 
characteristic unit for 

Equations (1-3) represent a vorticity-stream 
function formulation for the unsteady, 
incompressible Navier-Stokes equations. For 
numerical implementation, their component form in a 

1 2 

generalized coordinate system ( £ ,E ) has been 
given by Osswald, K. Ghia and U. Ghia [10]. Use of 
the two-dimensional component form of the equations 
1 2 

for oi and ii» in (E , E ) coordinates, together with 
appropriate boundary conditions for u and i|i, leads 
to a well-posed boundary-value problem for the 
airfoil flow. However the discussion on the 
boundary conditions will be deferred until the 
coordinate system to be used has been selected. 
The governing equations (1-3) form a coupled set of 
nonlinear equations and some comments are made next 
on the numerical method used to obtain their 
solution. 


3. Comments on the numerical Method 

The numerical method used was first developed 
by Osswald and Ghia [11] and was further refined by 
K. Ghia, Osswald and U. Ghia [1]. As discussed in 
Ref. [1], this method can be briefly described as 
follows. The spatial derivatives are approximated 
by appropriate finite-difference quotients using 
central differences for both convective and 
diffusive derivatives in the governing equations. 
It is significant to note that, In this study, even 
with central-diff erence approximations for all 
spatial derivatives, no artificial dissipation is 
added to dampen any high-frequency errors, but 
continued effort is made to carefully annihilate 
these errors through appropriate resolution of the 
various length scales of the problem. For 
consistent differencing, some of the metric 
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coefficients are evaluated at the staggered half- 
grid point locations, whereas the metric 
determinant, the solution field functions and the 
source term are evaluated at the cell corners. The 
conservative form of the vortlcity transport 
equation is solved, using an ADI method, whereas 
the stream function equation is solved using a 
block Gaussian elimination (BGE) technique. The 
BCE technique is a direct extension of the Gaussian 
elimination method to matrices whose individual 
elements are themselves matrices or blocks. The 
BGE technique provides the effective inversion of 
an (N-M x n-m) matrix through the actual inversion 
of a predetermined sequence of N(M«M) submatrices; 
the choice of M<N leads to the best computational 
efficiency. The BGE approach naturally divides 
itself into two separate calculation phases. In 
the first phase, a sequence of N(M«M) matrices is 
formed and individually inverted by simple scalar 
Gaussian elimination. This phase is the most time- 
consuming part of the calculation. Fortunately, 
this preliminary phase need be executed only once 
for a given coordinate choice, its result being 
permanently stored as a series of coefficient 
matrices for later use in the second phase of this 
procedure. The second phase consists of the actual 
solution of the block matrix problem for a 
prescribed source term through a set of recursion 
relationships. In this method, the transport 
equation at time level t +1 is discretized with the 

stream function being frozen at the time level t n - 

Due to this uncoupling, the formal temporal accuray 

of the scheme i3 0(At ) with the overall formal 
n 

truncation accuracy of 0[At n> (A? 1 ) 2 , (AC 2 ) 2 ]. For 

solution of the unsteady Navler-Stokes equations, 
it is the Dirichlet stream-function Poisson problem 
which must be solved repeatedly in a given 
coordinate system, for a progression of updated 
sources. Hence, the combined ADI-BGE method used 
here is very well suited for analyzing 
incompressible unsteady separated flows governed by 
the unsteady Navier-Stokes equations. 


4. Massive Separated Incompressible Flow 
Pa3t Symmetric Joukow3kl Airfoils 

The flow past a 12 percent thick symmetric 
Joukowski airfoil at moderate Re separates when the 
free stream is at incidence to the airfoil. As the 
angle of attack is increased, even at moderate Re, 
a region of massive separation develops on the 
upper surface. The geometry of the model problem 
of flow past a Joukowski airfoil can be accurately 
represented using conformal transformation 
techniques (see, e.g., Davis [12] and Ives [13]). 
and various conformal grids have been obtained for 
this problem. Furthermore, the flow past a 
Joukowski airfoil has been used by many 
investigators as a model problem for studying 
viscous separated flow, because of the simplicity 
of Its geometry for this class of flows. 
Consequently, results of other investigators are 
also available for qualitative comparison and, 
hence, this problem is used as a model problem in 
the present study. 

4 . 1 The Coordinate System 

Clustered conformal coordinates are employed 
in the present viscous-flow study. These 
coordinates are obtained by a parabolic 


transformation of the inviscid-f low complex- 
potential plane, followed by appropriate 1 -D 
clustering transformations for resolving the 
various length scales of the flow problem. 

For lifting cases, i.e., for flows past 
airfoils at non-zero incidence a, the corresponding 
inviscid flow has non-zero circulation T and, 
consequently, a multi-valued potential function 
$ in> Also, the inviscid flow is not symmetric 

about the airfoil mean line or camber line. 
Nevertheless, all of the viscous-flow results 
presented in this paper for symmetric, airfoils were 
obtained using a coordinate system which is 
symmetrical about the airfoil chord line, even for 
the lifting cases. In other words, the coordinates 
used are based on the corresponding potential flow 
with zero circulation and are obtained as a 
limiting case of the general coordinates developed 
for the inviscid flow with non-zero circulation. 
Some details of generating these general conformal 
coordinates for flows with circulation are 
described here in this section. For additional 
details, see Osswald, K. Ghia and U. Ghla [14]. 
The use of these general coordinates ,ln the study 
of viscous flows for lifting airfoils will be the 
subject of a future paper. 

4.1.1 Inviscid Flow Past a Joukowski Airfoil at 
Incidence 

The complex potential P for the inviscid flow 
of a uniform 3tream at Incidence o past a symmetric 
Joukowski airfoil is obtained conveniently from the 
knowledge of the complex potential for flow past a 
circular cylinder into which the Joukowski airfoil 
can be mapped via the Joukowski transformation. 
Accordingly, the complex potential is given as 

p ■ R + -r * 1 h in( !> (1,) 

where 

P » * 1 13 the complex potential, 

'^in’H'ln are the velocity potential and stream 

function, respectively, for the 2-D 
inviscid flow, 

A 1 2 

R - (r + lr ) is the complex variable in the 
circle plane in which the Joukowski 
airfoil profile is mapped as a 
circle, 

a Is the radius of the circular cylinder in the 
circle plane, 
and 

r • 4ira sina is the circulation around the 
cylinder (or airfoil). 

The radius a is related to the thickness ratio 

(i.e. maximum thickness t /chord c) of the 

max 

airfoil and the incidence angle a of the free 
stream approaching the airfoil. This relation is 
established by the Joukowski transformation as 
discussed next. 

The complex variable R in the circle plane is 

related to the complex variable z (= x^ix 2 ) in the 
physical plane via the Joukowski transformation as 
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(5) 


where 


o la 

Re - ep , 


LE 1 


C<T2e> ♦ C^)] p . 


p - ( 1+2e)/[1 + (1+2e)(3+2e)] 

and e Is related to the thickness ratio t^^/c by 
the expression 


4ep( 1 +e ) 2 

1 + 2e( 1 +e ) ( 1 -cose ) 
ro 


t 


[sine (1-eos8 )] . 

c m m 


( 6 ) 


In this expression, 9 is the value of 8 in the 
m 

physical plane at which the maximum thickness of 
the airfoil occurs, i.e., 


In fact, for each point on the trailing streamline 
in the physical plane. 



This distinction in the values of 4, must be 

l n 

maintained for all points along the upper and lower 
trailing streamlines, downstream of the trailing 
edge. Hence, a desired grid-point distribution may 
be assigned only along the upper trailing 
streamline, for example, while the grid point 
distribution along the lower trailing streamline is 
no longer arbitrary, but is constrained by 
condition (8). 

Each of the above three points 13 addressed as 
follows. First of all, the P-plane i3 transformed 

to yet another plane, namely, the n-plane (n - n 1 
2 

+in )using the conformal parabolic transformation 


9 m ' 6 }y<S)-0-5 t MX /c • 


2 

It should be noted that' an 0(e ) approximate 
solution to Eq. (6) can be obtained to yield 

8 = 2ir/3 so that e => — — -~ x . In the present 

ra 3/3 ° 

study, however, this approximation was not invoked, 
but Eq. (6) was solved numerically using the method 
of successive substitution. It was observed that, 
for a 12 percent thick symmetric Joukowski airfoil, 
the value of e so determined differed from the 
approximate value by almost 10 percent. 


Figure (1) shows the inviscid flow problem 

considered in the physical plane, the circle plane 

and, finally, In the complex-potential plane. The 

lines of 4. » constant and ¥. » constant form a 

in in 

rectangular coordinate system in the P-plane and 
yield a surface-oriented conformal coordinate 
system in the physical plane. Hence, the P-plane 
could constitute the computational plane. However, 
it is not quite appropriate for this purpose 
because of the following three main reasons. First 

of all, the coordinates ¥, » constant experience 

in 

large turning in the vicinity of the leading 
stagnation point, the degree of turning depending 
upon the incidence angle a and the roundedness of 
the airfoil in this region and near the leading 
edge. Also, the manner in which the coordinates 
are disposed is not conducive to proper resolution 
in this region. Secondly, the flow occurs over the 
entire unbounded P-plane extending from -« to +» 

along both 4, and ¥, ; this is computationally 
in in 

impracttal. Finally, it Is noted that the airfoil 
profile between the leading stagnation point (LSP) 
and the trailing edge point (TEU) following the 


upper contour, i.e., along , ? ln ■ 0 + , extends 

from 4. = 0 to 4, = 4_„.. Following the lower 

in m i t.u 


contour, i.e., 
between LSP and 
4 tel where 

*TEL “ ^TEU “ 


along ¥ in = 0 , the airfoil profile 
TEL extends from 4 - 0 to 4, » 


r . 


(7) 


2 

(n) - P + 2a[cosa + (a+ir)sina] . (9) 

In the n-plane, the flow problem occupies only the 
upper-half plane. Secondly, suitable 1-D 
clustering transformations are employed to map the 
upper half of the n-plane to a unit square In the 

5-plane (5 » 5 1 ♦ i£ 2 ) while also providing 
resolution of the various streamwise and normal 
length scales in the flow. The nature of these 
clustering transformations is as discussed in 
Ref. [15]. Differences In details are necessitated 
by facts such as the leading stagnation point LSP 
no longer coinciding with the leading edge point 
LE. Finally, the requirement expressed in Eqs. (7) 
or (8) is implemented in the n-plane by recognizing 
that, downstream of the trailing edge point, 

along * ln » 0 + , P = 4 u , 

and, along t • 0 , P = 4^ , 
so that 

(n l )2 " ( \ )2 " r HO) 

for points downstream of the trailing edge. This 
implementation is carried over to the 5-plane as 
well. The 5-plane constitutes the final 
computational plane. Use of a uniform mesh in the 
unit square in the 5-plane provides a surface- 
oriented mesh In the unbounded physical plane, with 
appropriately clustered grid-point distributions. 
Figure 1(d-e) shows the flow problem represented in 
the n-plane and finally, in the 5~plane. 

The procedure described in this section was 
used to generate clustered conformal coordinates 
for analysis of viscous flow past a 12 percent 

thick (i.e., t /c » 0.12) Joukowski airfoil, 
max 

Figure 2(a) shows the coordinate system 
corresponding to a (229 x 45) mesh for a free 
stream at zero incidence, i.e., a ° 0, while Fig. 
2(b) shows the coordinates for the case with 
a = 10°. Clustering of coordinate lines can be 
observed in the region near the airfoil surface as 
well as near the leading and trailing edges and the 
leading stagnation point. 
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4.2 Initial and Boundary Conditions 

The final working variables in the flow- 
problem are the vortlcity u and the perturbation 
stream function i|> defined as the departure of the 

actual viscous-flow stream function from the 
corresponding invlscid solution ¥ , i.e. 


The initial conditions used for the unsteady 
solution correspond to the requirement that 
invlscid flow prevails everywhere at the starting 
time instant T » 0. 

The boundary conditions u sed correspond to the 
condition of zero slip at the airfoil surface, 
inviscld-flow conditions infinitely far from the 
airfoil and continuity of the flow solution across 
the branch cut along the trailing streamline in the 
wake. The actual numerical form of the second- 
order vortlcity boundary condition used at the 
surface is given in Ref. [1]. As in the previous 
studies by the authors, implicit treatment of all 
boundary conditions is emphasized. This is 
maintained for the branch-cut boundary as well. 


5. Results and Discussion 
5.1 Physical Scales 

The unsteady Hav ier-Stokes analysis and 
solution procedure discussed in the earlier 
sections are applied to a 12 percent thick 
symmetric Joukowski airfoil at moderate Re and 
angle of incidence a. The steady separated or the 
unsteady massive separated flow fields over the 
Joukowski airfoil at angle of incidence have 
multiple disparate length and time scales. If the 
physics of this class of flow problems has to be 
accurately modelled, then these scales must be 
resolved; otherwise, the flow field will not be 
simulated accurately in critical areas such as near 
the separation and reattachment points, etc. 
Rothmayer and Davis [16] have discussed how the 
lack of proper resolution of the length scale at 
the separation point could rather drastically 
influence the overall eddy and the reattachment 
point if there is one. Based on whether the 
separation is mild or massive, two distinct flow 
structures result. For the former, the viscous 

layer is of 0(t/c), where t i3 the thickness of the 
airfoil, whereas, for the latter case, the detached 
zone is of 0(1). Some of the important length 
scales associated with these flow structures are 

- 1/2 

the boundary- layer scale of 0(Re ), the 

separation scales of 0(Re J ) and 0(Re J ) in the 
streamwise and normal directions, respectively, at 
the separation point. Further, since the 
boundaries external to the body are at true 
infinity, the invisoid scale of 0(u) must be 
resolved so that the flow smoothly asymptotes to 
the invlscid boundary conditions. In order to 
follow the vortices in the wake region, the length 
scale of the wake region must be resolved. For the 
unsteady flow field, the temporal scale based on 
the highest critical frequency needs to be 
resolved. In the scope of the present study, the 
scale based on the frequency of the self-induced 
unsteadiness generated by the shedding of vortices 
must be resolved. Only if all of the scales 


mentioned are properly considered in the generation 
of the grid is there a hope of simulating physical 
phenomena using the unsteady Navier-Stokes 
equations. 

5.2 Quality of Grid, Accuracy and Computational 

Efficiency 

This discussion on the physical scales shows 
large disparity in the various scales. This places 
a major burden on the grid. The conformal 
clustered grid shown in Fig. 2(a,b) has different 
clustering transformations on the upper and lower 
surfaces of the airfoil in the streamwise direction 
and a different clustering in the wake. Two 
distinct normal transformations are used both on 
the upper and lower surfaces. The transformation 
metrics are smooth and their continuity is 
maintained across the various regions. Even after 
exercising. this care, a (229x45) computational grid 
was considered necessary; it uses up the, full core 
capacity of the present host AMDAHL 470 V/7A and 
Perkin Elmer 3250 MPS computer systems. Hence, the 
standard approach of refining the grid is not 
feasible since computational resources do not 
permit it. Therefore, the results of thi3 study, 
at present, are considered only qualitative. Every 
effort is made to prevent them from being 
contaminated by numerical dissipation or quality of 
grid; but only after thorough comparison with 
experimental or asymptotic studies can their 
quantitative usefulness be assessed. Even in the 
present work, the analysis was checked for 
consistancy by using the free-stream condition at 
all boundaries and ensuring the solution recovered 
the free-stream values everywhere. It is 
anticipated that the U3e of the general conformal 
clustered grid, discussed in section 4.1, for which 
the wake centerline approaches the free stream 
direction will enhance the accuracy of the results. 

The use of central differences throughout the 
flow field leads to the formal overall accuracy of 

0[ At , ( A6 1 ) 2 ,(A£ 2 ) 2 ]. For the configurations with 
Re = 1000, the resulting solutions are wiggle-free; 
however, for the case of Re - 10,000, the 
inadequacy of the grid leads to some oscillations 
In the vortlcity; these are presently being 
examined. The relative computational efficiency of 
the present algorithm was measured in terms of the 
CPU time t required to advance the solution by one 
time step per spatial grid point. The value of the 

-4 

"computational effort" index is t « 3-72 x 10 
seconds for the AMDAHL 470 V/7A computer. 

5.3 Mildly Separated Steady Flow: Re -1000, 



The flow past a symmetric Joukowski airfoil 
with the parameters listed was computed starting 
from the Invlscid solution shown in Fig. 3(a) at 
the characteristic time of T = 0. This, as well 
as, the other results to follow, have been computed 
using a coordinate system which is symmetrical 
about the airfoil chord line and is shown in Figure 
3(b). The transient stream function and vortlcity 
contours leading to steady state are shown in Fig. 
3(0-1) . At T = 1.92, in Fig. 3(e), a 3mall 
separation bubble has been established near the 
trailing edge (TE) . This primary bubble bursts, 
and the shear layer detaches, with a simultaneous 
emergence of a secondary bubble at the TE, as seen 
in Fig. 3(g). The secondary bubble intensifies and 
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grows and merges with the primary bubble as seen In 
Fig. 3(1),. This configuration achieves steady 
state around T = 7, but the calculations were 
continued up to T = 13, Fig. 3(k), by which time It 
became certain that steady state has been achieved. 
The vortlclty contours show a maximum vorticity of 
296 at T * 0.4, Fig. 3(d); the maximum vorticity 
gradually reduces to 222 at T > 13.0, Fig. 3(1). 
The vorticity contours show laminar boundary layers 
on the upper and lower surfaces with a tongue-like 
behavior in the vorticity at T ■ 1.92 in Fig. 3(f) 
and at subsequent times. As time progresses, an 
asymmetry developes in the vorticity contours and 
it appears that the region of concentrated 
vorticity near the leading edge (LE) has been 
satisfactorily resolved. 

5. 1 * Massive Separa tion - Unstead y Flow; 

Re - 1000, a =» 15° 

Although the present results are obtained as 
the solution of the unsteady Mav i er -S t o k e s 
equations, the point of flow separation has been 
defined as the point on the airfoil surface where 
the surface shear vanishes and the instantaneous 
streamline detaches from the solid wall. (For 
further details regarding the separation point in 
unsteady flow, the reader is referred to the review 
article by McCrowskey [17]). For the flow 
configuration discussed in this subsection, results 
from T - 10.0 to T ■ 21, extending up to two chord 
lengths in the wake region, are shown in Fig. 4(a- 
f). Results beyond the characteristic time of 
T » 21 , not shown here, confirm that the flow is 
inherently unsteady due to the vortex shedding 
mechanism established by the free shear-layer 
instability. The instantaneous stream function 
contours at T = 10 show massive separation on the 
suction surface. The shear layer detaches near the 
LE and never reattaches to the airfoil surface. 
Two counterclockwise co-rotating bubbles appear aft 
of the shoulder towards the TE. These bubbles are 
in the process of coalescing. The corresponding 
instantaneous vorticity contours show several well 
developed eddies, some of which are in the process 
of merging due to intense interaction between them. 
Intermediate solutions are needed to explain the 
evolution process of the various eddies and the 
collective Interaction phenomenon described by 
Ho et al. [18], where more than two eddies merge to 
form a single eddy. The flow field is dominated by 
large vortices and their passage over the upper 
surface produces unsteady forces that are different 
from those observed during static stall. Following 
Smith [9], Rothmayer and Davis [16] have referred 
to this phenomenon as 'dynamic stall'. The large- 
scale organized structure observed needs further 
examination in order to seek its relation, if any, 
to turbulent coherent structure (see Hussain [19]). 
Figure 4(c) shows the emergence of a new bubble on 
the suction surface and the previous bubbles in 
Fig. 4(a) appear to have merged and convected 
downstream. The corresponding vorticity contours 
at T » 11 are shown in Fig. 4(d) and are 
characterized by very similar vortical structures 
as in Fig. 4(b). At T • 21 , the flow structure in 
Fig. 4(e-f) resembles that in Fig. 4(c-d). 

For this configuration, the results are also 
plotted over a wider field extending up to five 
chord lengths in the wake region and are shown in 
Fig. 5(a-d) at the characteristic times T = 11 and 
21. The massive separated unsteady flow field is 
clearly seen here in Fig. 5(a-c). In Fig. 5(c), 
the shear layer that detaches near LE reattaches at 


the TE. It is Important that instantaneous 
coefficients of lift and drag be calculated and 
correlaton with streamline contours be established 
to show the increase in lift coefficient that may 
result due to the reattachment of the shear layer 
at the TE. Figure 5(b-d) shows several vortices 
interacting with each other. Some fragmentation of 
lage-scale eddies into small eddies can also be 
observed. Further analysis of the results for this 
case reveals a nonlinear limit-cycle type behavior 
which is presently being quantified, so as to 
determine, for instance, the Strouhal number for 
vortex shedding process. 

The velocity vectors and pressure coefficient 
0^ for the two Joukowski airfoil configurations 

studie at Re = 1000 are depicted in Fig. 6(a-d). 
For the case of a = 5°, steady-state velocity 
profiles at various locations on the airfoil and in 
the wake are shown in Fig. 6(a). These profiles 
clearly show the mildly separated region as well as 
the velocity deficit in the wake region. The 
inviscid as well as the viscous C p are shown in Fig. 

6 (b); the maximum difference between them occurs 
on the suction surface, very close to the LE, with 
the magnitude of the inviscid C p being considerably 

larger than that of the viscous C p , as expected. 

For the unsteady massive separated flow case, 

instantaneous velocity vectors as well as C are 

P 

shown in Fig. 6(c-d) at T « 21. The velocity 
profiles show the massive separated region on the 
airfoil surface; the presence of multiple bubbles 
in the wake can also be 3een by the nature of the 
wake velocity profiles. The curves of 

instantaneous C show that, for the viscous flow, 
P 

the leading stagnation point has shifted towards 
the LE. 

The computation of this unsteady flow case 
with massive separation demonstrates the ability of 
the analysis to treat the flow past this Joukowski 
airfoil at higher incidence. Elrafee, Wu and 
Lekoudis [20] have used a 9 percent thick Joukowski 
airfoil at Re - 1000, a - 15° and have computed 
subsonic flow, at M • 0.4 and Pr « 1.0, around it. 
The results for instantaneous streamline and 
vorticity contours, as well as the C p distribution, 

show rather minimal departure from the 
corresponding incompressible flow. Further, 
Sugavanam and Wu [21] have computed turbulent flow 
past a modified 12 percent thick Joukowski airfoil 

at Re « 3.6*10^, a => 15°, using a two-equation k-e 
model. The contours of time-averaged streamlines 
show qualitative resemblance to the contours of 
instantaneous streamlines shown in Fig. 4, except 
that the separated region at the high-Re case is 
smaller in extent. It would be interesting to 
qualitatively compare the time-averaged vorticity 
contours with those in Fig. 4; however, these have 
not been presented by those authors. It should be 
noted that, in their work, the far-field boundary 
conditon was placed at approximately eight chord 
lengths. Detailed examination of other parameters 
which are more sensitive could perhaps better 
reveal the departure from the present 
Incompressible unsteady case computed. 
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5.5 M as sive Separation -Unsteady Flow: 

Re -■ IoTQOO. a - 5° 

McCroskey [17] has discussed that, when 
Re > 1000, three-dimensional and turbulence effects 
are present in the flow field and the unsteady 
Navier-Stokes analysis should account for these 
effects. In the present study, this configuration 
is used to test the ability of the code to compute 
this flow as well as to find the features 
distinguishing this flow from that for Re = 1000, 
even if they are only qualitative, to aid in future 
study of this phenomenon. As stated earlier, the 
symmetric grid used is inadequate in the normal 
direction as well as in the wake, since the wake 
centerline does not follow the coordinates used. 
Also, the vorticity contours display some wiggles. 
Hence, the results for this configuration are, at 
best, qualitative. Figure 7(a-£) 3hows the 
Instantaneous stream function and vorticity 
contours from T « 2.0 to T - 34 for this unsteady 
flow. At T ■ 2.0, the instantaneous streamlines in 
Fig. 7(a) show a well-behaved laminar separated 
flow. From the corresponding vorticity contours in 
Fig. 7(b), it appears that the boundary layer is 
well behaved and a narrow wake trails the airfoil. 
By T - 5.0, Fig. 7(c), most of the suction surface 
exhibits reversed flow. Figure 7(d) shows the 
onset of Tollmien-Schlichting type Instability on 
the upper surface, as well as the emergence and 
intensification of an eddy near the trailing edge. 
At present, the available graphics facility is not 
adequate for post-processing a large-scale data 
base and only painstaking hard labor has permitted 
generation of computer plots at some characteristic 
time instants. With proper graphics facilities 
which are presently being sought, more precise 
informaton will be made available in the future. 
Somewhere between T * 5.0 and T - 12.0, a total 
breakdown of the flow occurs, with a sudden 
Increase in the normal length scale as can be seen 
on the suction surface in Fig. 7(f). The 
streamline contours in the wake at T - 12.0 in Fig. 
7(e) are not smooth; this is due to the wiggles 
that appear in the vorticity field. The wiggles 
are more prominent on the wake originating from the 
lower surface and are somewhat reduced as T 
increases to 14.0, Fig. 7(h), and to subsequent 
times at T * 24.0 and T - 34, as depicted in Fig3. 
7(i) and 7(k). The streamline contours at T * 14.0 
in Fig. 7(i) bear a strong similarity to those at 
T = 24.0, suggesting that the flow may exhibit a 
nonlinear limit-cycle behavior. These qualitative 
results are encouraging and warrant careful 
investigation of high-Re flows using this unsteady 
analysis. 


6. Conclusions 

The unsteady analysis of Osswald and Ghia [11] 
has been extended to analyze 2-D unsteady separated 
external flow past symmetric Joukowski airfoils at 
high incidence and moderate Re. The boundaries 
external to the airfoils are placed at infinity. 
The discretized problem is formulated using central 
differences for spatial derivatives, thus avoiding 
any artificial viscosity. The fully implicit ADI- 
BGE time-marching method, with formal overall 

accuracy of 0[&t, (AC 1 ) 2 , (AC 2 ) 2 ], is used to solve 
the discretized equations. 

Three configurations are investigated for a 12 
percent thick Joukowski airfoil and their flow 


features are carefully discussed. For Re - 1000 
and a =5°, a steady separated flow structure is 
obtained. For this Re, when a is increased to 15°, 
a massive unsteady separated flow field is 
obtained. A nonlinear limit-cycle type analysis is 
attempted on the latter results. At angle of 
incidence of 5°, Re was increased to 10,000 and the 
flow shows an Instability around T = 5.0, and 
exhibits a turbulent-like behavior thereafter. The 
results of the present analysis are still 
considered qualitative and fine-mesh results using 
the new clustered conformal grid (Fig. 2b), are 
desired in order to make more conclusive statements 
about the flow structure observed in these results. 
The results obtained do show the potential of the 
present analysis to study high-incidence high-Re 
flow. It is planned to extend the analysis to 
lifting NACA airfoils and carefully compare the 
results with available experimental and numerical 
results. 
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(e) Comoutational Plane 
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Fig. 1. Representation of Inviscid Flow Past Symmetric JouXowski Airfoil with Circulation, 
in Physical Plane and Various Transformed Planes. 
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INSTANTANEOUS STREAM-FUNCTION CONTOURS INSTANTANEOUS VORTICITY CONTOURS 

Fig. 3. Flow Past Symmetric Joukowski Airfoil at Re = 1,000, a=5°. T=0.40 and 1 . 12 . 

[At T=0: Inviscid Stream Function and Grid-Point Distribution) 
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Fig. 5. Flow Past Symmetric Joukowski Airfoil at Re = 1,000, a = 15°. T = 11.00 and 21.00 
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Velocity Vectors and Surface Pressure Coefficient for Flow Past Symmetric Joukowski Airfoil 
at Incidence; Re = 1,000. 



iscous STRenn function 



NSTANTANEOUS STREAM-FUNCTION CONTOURS INSTANTANEOUS VORTICITY CONTOURS 








